igl, MeshCAD, Polyfem 3D example
import numpy as np import math import igl import meshplot as mp import wildmeshing as wm import polyfempy as pf
Read the mesh
V, F = igl.read_triangle_mesh("mesh.obj") mp.plot(V, F)
Get mesh bounding box and center
minn = np.min(V, axis=0) maxx = np.max(V, axis=0) center = (minn+maxx)/2
Function to select sidesets:
- Sideset 2: z-component close to the bottom boundary
minn[2]and larger than y of center - Sideset 3: z-component close to the bottom boundary
minn[2]and smaller than y of center
def sideset(p): if abs(p[2] - minn[2]) < 0.5: if p[1] > center[1]: return 2 else: return 3 return 1
Wildmeshing tetrahedralization
wm.tetrahedralize("mesh.obj", "out.mesh", mute_log=True)
Loading the mesh
solver = pf.Solver() solver.load_mesh_from_path("out.mesh")
[2019-07-22 19:36:29.077] [polyfem] [info] Loading mesh... [2019-07-22 19:36:29.078] [geogram] [info] Loading file out.mesh... [2019-07-22 19:36:29.140] [geogram] [info] (FP64) nb_v:4019 nb_e:0 nb_f:4872 nb_b:0 tri:1 dim:3 [2019-07-22 19:36:29.140] [geogram] [info] nb_tets:15744 [2019-07-22 19:36:29.140] [geogram] [info] Attributes on vertices: point[3] [2019-07-22 19:36:29.326] [polyfem] [info] mesh bb min [-92.075, -98.5272, -42.037], max [92.075, 28.6776, 42.037] [2019-07-22 19:36:29.327] [polyfem] [info] took 0.249492s
Using the previous function as sidesets and visualization for confirmation
solver.set_boundary_side_set_from_bary(sideset) ps, ts, s = solver.get_boundary_sidesets() col = np.zeros_like(s) col[s==2] = 2 col[s==3] = 3 mp.plot(ps, ts, col, shading={"wireframe": False})
Setup the problem
settings = pf.Settings() problem = pf.Problem() settings.set_pde(pf.PDEs.LinearElasticity) settings.set_material_params("E", 200) settings.set_material_params("nu", 0.35) problem.set_displacement(2, [0, 0, 0]) problem.set_force(3, [0, 0.5, 0]) settings.set_problem(problem)
Solve!
solver.settings(settings) solver.solve()
[2019-07-22 19:36:29.474] [polyfem] [info] simplex_count: 15744 [2019-07-22 19:36:29.474] [polyfem] [info] regular_count: 0 [2019-07-22 19:36:29.474] [polyfem] [info] regular_boundary_count: 0 [2019-07-22 19:36:29.474] [polyfem] [info] simple_singular_count: 0 [2019-07-22 19:36:29.474] [polyfem] [info] multi_singular_count: 0 [2019-07-22 19:36:29.474] [polyfem] [info] boundary_count: 0 [2019-07-22 19:36:29.474] [polyfem] [info] multi_singular_boundary_count: 0 [2019-07-22 19:36:29.474] [polyfem] [info] non_regular_count: 0 [2019-07-22 19:36:29.474] [polyfem] [info] non_regular_boundary_count: 0 [2019-07-22 19:36:29.474] [polyfem] [info] undefined_count: 0 [2019-07-22 19:36:29.474] [polyfem] [info] total count: 15744 [2019-07-22 19:36:29.475] [polyfem] [info] Building isoparametric basis... [2019-07-22 19:36:29.528] [polyfem] [info] Computing polygonal basis... [2019-07-22 19:36:29.528] [polyfem] [info] took 0.000912479s [2019-07-22 19:36:29.531] [polyfem] [info] hmin: 0.741093 [2019-07-22 19:36:29.531] [polyfem] [info] hmax: 18.3008 [2019-07-22 19:36:29.531] [polyfem] [info] havg: 5.9001 [2019-07-22 19:36:29.531] [polyfem] [info] took 0.0526661s [2019-07-22 19:36:29.531] [polyfem] [info] flipped elements 0 [2019-07-22 19:36:29.531] [polyfem] [info] h: 18.3008 [2019-07-22 19:36:29.531] [polyfem] [info] n bases: 4019 [2019-07-22 19:36:29.531] [polyfem] [info] n pressure bases: 0 [2019-07-22 19:36:29.531] [polyfem] [info] Assigning rhs... [2019-07-22 19:36:29.541] [polyfem] [info] took 0.0102129s [2019-07-22 19:36:29.541] [polyfem] [info] Assembling stiffness mat... [2019-07-22 19:36:29.674] [polyfem] [info] took 0.13286s [2019-07-22 19:36:29.674] [polyfem] [info] sparsity: 435691/145371249 [2019-07-22 19:36:29.674] [polyfem] [info] Solving LinearElasticity with [2019-07-22 19:36:29.676] [polyfem] [info] Hypre...
Visualization of displacement and stresses
p, t, d = solver.get_sampled_solution() m, _ = solver.get_sampled_mises_avg()
mp.plot(p+d, t, m, shading={"wireframe": False})